Benchmarking machine learning and parametric methods for genomic prediction of feed efficiency-related traits in Nellore cattle

Genomic selection (GS) offers a promising opportunity for selecting more efficient animals to use consumed energy for maintenance and growth functions, impacting profitability and environmental sustainability. Here, we compared the prediction accuracy of multi-layer neural network (MLNN) and support vector regression (SVR) against single-trait (STGBLUP), multi-trait genomic best linear unbiased prediction (MTGBLUP), and Bayesian regression (BayesA, BayesB, BayesC, BRR, and BLasso) for feed efficiency (FE) traits. FE-related traits were measured in 1156 Nellore cattle from an experimental breeding program genotyped for ~ 300 K markers after quality control. Prediction accuracy (Acc) was evaluated using a forward validation splitting the dataset based on birth year, considering the phenotypes adjusted for the fixed effects and covariates as pseudo-phenotypes. The MLNN and SVR approaches were trained by randomly splitting the training population into fivefold to select the best hyperparameters. The results show that the machine learning methods (MLNN and SVR) and MTGBLUP outperformed STGBLUP and the Bayesian regression approaches, increasing the Acc by approximately 8.9%, 14.6%, and 13.7% using MLNN, SVR, and MTGBLUP, respectively. Acc for SVR and MTGBLUP were slightly different, ranging from 0.62 to 0.69 and 0.62 to 0.68, respectively, with empirically unbiased for both models (0.97 and 1.09). Our results indicated that SVR and MTGBLUBP approaches were more accurate in predicting FE-related traits than Bayesian regression and STGBLUP and seemed competitive for GS of complex phenotypes with various degrees of inheritance.


Material and methods
The FE-related traits and genomic information were obtained for 1,156 animals from an experimental breeding program at the Beef Cattle Research Center (Institute of Animal Science -IZ).

Phenotypic and genotypic information
Animals were from an experimental breeding program at the Beef Cattle Research Center at the Institute of Animal Science (IZ) in Sertãozinho, São Paulo, Brazil.Since the 1980s, the experimental station has maintained three selection herds: Nellore control (NeC) with animals selected for yearling body weight (YBW) with a selection differential close to zero, within birth year and herd, while Nellore Selection (NeS) and Nellore Traditional (NeT) animals are selected for the YBW with a maximum selection differential, also within birth year and herd 25 .In the NeT herd, sires from commercial herds or NeS eventually were used in the breeding season, while the NeC and NeS were closed herds (only sires from the same herd were used in the breeding season), with controlled inbreeding rate by planned matings.In addition, the NeT herd has been selected for lower residual feed intake (RFI) since 2013.In the three herds, the animal selection is based on YBW measured at 378 days of age in young bulls.
The FE-related traits were evaluated on 1156 animals born between 2004 and 2015 in a feeding efficiency trial, in which they were either housed in individual pens (683 animals) or group pens equipped with the GrowSafe feeding system (473 animals), with animals grouped by sex.From those, 146 animals were from the NeC herd (104 young bulls and 42 heifers), 300 from the NeS herd (214 young bulls and 86 heifers), and 710 from the NeT herd (483 young bulls and 227 heifers).Both feeding trials comprised at least 21 days for adaptation to the feedlot diet and management and at least 56 days for the data collection period.The young bull and heifers showed an average age at the end of the feeding trial was 366 ± 27.5 and 384 ± 45.4 days, respectively.
A total of 780 animals were genotyped with the Illumina BovineHD BeadChip assay (770 k, Illumina Inc., San Diego, CA, USA), while 376 animals were genotyped with the GeneSeek Genomic Profiler (GGP Indicus HD, 77 K).The animals genotyped with the GGP chip were imputed to the HD panel using FImpute v.3 26 with an expected accuracy higher than 0.97.Autosomal SNP markers with a minor allele frequency (MAF) lower than 0.10 and a significant deviation from Hardy-Weinberg equilibrium (P ≤ 10 −5 ) were removed, and markers and samples with call rate lower than 0.95 were also removed.An MAF lower than 10% was used to remove genetic markers with lower significance and noise information in a stratified population.After this quality control procedure, genotypes from 1,024 animals and 305,128 SNP markers remained for GS analyses.Population substructure was evaluated using a principal component analysis (PCA) based on the genomic relationship matrix using the ade4 R package (Supplementary Figure S1) 27 .

Feed efficiency-related traits
Animals were weighed without fasting at the beginning and end of the feeding trial, as well as every 14 days during the experimental period.The mixed ration (dry corn grain, corn silage, soybean, urea, and mineral salt) was offered ad libitum and formulated with 67% of total digestible nutrients (TDN) and 13% of crude protein (CP), aiming for an average daily gain (ADG) of 1.1 kg.
The following feed efficiency-related traits were evaluated: ADG, dry matter intake (DMI), feed efficiency (FE), and RFI.In the individual pens, the orts were weighed daily in the morning before the feed delivery to calculate the daily dietary intake.In the group pens, the GrowSafe feeding system automatically recorded the feed intake.Thus, the DMI (expressed as kg/day) was estimated as the feed intake by each animal with subsequent adjustments for dry matter content.ADG was estimated as the slope of the linear regression of body weight (BW) on feeding trial days, and the FE was expressed as the ratio of ADG and DMI.Finally, RFI was calculated within each contemporary group (CG), as the difference between the observed and expected feed intake considering the average metabolic body weight (MBW) and ADG of each animal (Koch et al., 1963) as follows: where β 0 is the model intercept, β 1 and β 2 are the linear regression coefficients for ADG and MBW = BW 0.75 , respectively, and ε is the residual of the equation representing the RFI estimate.
The contemporary groups (CG) were defined by sex, year of birth, type of feed trial pen (individual or collective) and selection herd.Phenotypic observations with values outside the interval of ± 3.5 standard deviations below and above the mean of each CG for each trait were excluded, and the number of animals per CG ranged from 10 to 70.

Estimation of genetic parameters
The (co)variance components and heritability for FE-related traits were estimated considering a multi-trait GBLUP (MTGBLUP) as follows: Where y is the matrix of phenotypic FE-related traits (ADG, FE, DMI, and RFI) of dimension Nx4 ( N individuals and four traits); β is the vector of fixed effects, linear and quadratic effects of cow age, and linear effect of animal's age at the beginning of the test; a is the vector of additive genetic effects (breeding values) of animal, and e is a vector with the residual terms.The X and Z are the incidence matrices related to fixed (b) and random effects (a), respectively.It was assumed that the random effects of animals and residuals were normally distributed, as a ∼ N(0, G ⊗ S a ) and e ∼ N(0, I ⊗ S e ), where G is the additive genomic relationship matrix between genotyped individuals according to VanRaden 28 , I is an identity matrix, ⊗ is the Kronecker product, and  are the additive genetic and residual (co)variance matrices, respectively.The G matrix was obtained according to VanRaden 28 : where M is the SNP marker matrix with codes 0, 1, and 2 for genotypes AA, AB, and BB adjusted for allele frequency expressed as 2p j , and p j is the frequency of the second allele jth SNP marker.The analyses were performed using the restricted maximum likelihood (REML) method through airemlf90 software 29 .The predictf90 software 29 was used to obtain the phenotypes adjusted for the fixed effects and covariates ( y * = y − X β ).The adjusted phenotypes were used as the response variable in the genomic predictions.
Tthe GEBVs accuracy ( Acc GEBV ) in the whole population, was calculated based on prediction error variance (PEV) and the genetic variance for each FE-related trait ( σ 2 a ) using the following equation 30 :

Data splitting
A forward validation scheme was applied for computing the prediction accuracies using machine learning and parametric methods, splitting the dataset based on year of birth, with animals born between 2004 and 2013 assigned as the reference population (n = 836) and those born in 2014 and 2015 (n = 188) as the validation set.
For ML approaches, we randomly split the training dataset into fivefold to train the models.

Parametric methods
Single-trait GBLUP (STGBLUP).Genomic prediction for FE-related traits considering the STGBLUP can be described as follows: where y * is the Nx1 vector of adjusted phenotypic values for FE-related traits, µ is the model intercept, Z is the incidence connecting observations; a is the vector of predicted values, assumed to follow a normal dis- tribution given by N(0, Gσ 2 a ) and e is the Nx1 vector of residual values considered normally distributed as N(0, Iσ 2 e ), in which I is an identity matrix, σ 2 e is the residual variance.The STGBLUP model was performed using blupf90 + software 29 .

Multi-trait GBLUP (MTGBLUP)
Genomic prediction for FE-related traits considering MTGBLUP can be described as follows: where y * is the matrix of adjusted phenotypes of dimension Nx4, µ is the trait-specific intercept vector, Z is the incidence matrix for the random effect; a is an Nx4 matrix of predicted values, assumed to follow a normal distribution given by MVN(0, G ⊗ S a ) where S a represents genetic (co)variance matrix for the FE-related traits (4 × 4).The residual effects (e) were considered normally distributed as MVN(0, I ⊗ S e ) in which I is an identity matrix, and S e is the residual (co)variance matrix for FE-related traits (4 × 4).The MTGBLUP was implemented in the BGLR R package 14 considering a Bayesian GBLUP with a multivariate Gaussian model with an unstructured (co)variance matrix between traits ( S a ) using Gibbs sampling with 200,000 iterations, including 20,000 samples as burn-in and thinning interval of 5 cycles.Convergence was checked by visual inspection of trace plots and distribution plots of the residual variance.

Bayesian regression models
Five Bayesian regression models with different priors were used for GS analyses: Bayesian ridge regression (BRR), Bayesian Lasso (BL), BayesA, BayesB, and BayesC.The Bayesian algorithms for GS were implemented using the R package BGLR version 1.09 14 .The BGLR default priors were used for all models, with 5 degrees of freedom (dfu), a scale parameter (S), and π.The Bayesian analyses were performed considering Gibbs sampling chains of 200,000 iterations, with the first 20,000 iterations excluded as burn-in and a sampling interval of 5 cycles.Convergence was checked by visual inspection of trace plots and distribution plots of the residual variance.For Bayesian regression methods, the general model can be described as follows: where µ is the model intercept; x iw is the genotype of the ith animal at locus w (coded as 0, 1, and 2); u w is the SNP marker effect (additive) of the w-th SNP (p = 305,128); and e i is the residual effect associated with the observation of ith animal, assumed to be normally distributed as e ∼ N(0, Iσ 2 e ).The BRR method 14 assumes a Gaussian prior distribution for the SNP markers ( u w ) , with a common variance (σ 2 u ) across markers so that p u 1 , . . .
and the residual variance is also assigned a scaled-inverse Chi-squared distribution with degrees of freedom (df e ) and scale parameters (S e ).
Bayesian Lasso (BL) regression 31 used an idea from Tibshirani 32 to connect the LASSO (least absolute shrinkage and selection operator) method with the Bayesian analysis.In the BL, the source of variation is split into residual term (σ 2 e ) and variation due to SNP markers ( σ 2 u w ).The prior distribution for the additive effect of the SNP marker p u w |τ 2 j , σ 2 e follows a Gaussian distribution with marker-specific prior variance given by p u w |τ 2 j , σ 2 e = p w=1 N(u w 0, τ 2 j σ 2 e .This prior distribution leads to marker-specific shrinkage of their effect, whose their extent depends on the variance parameters τ 2 j .The variance parameters τ 2 j is assigned as exponential independent and identically distributed prior, p τ Exp τ 2 j 2 and the square lambda regularization parameter ( 2 ) follows a Gamma distribution ( p 2 = Gamma(r, θ)) , where r and θ are the rate and shape parameters, respectively 31 .Thus, the marginal prior for SNP markers is given by a double exponential (DE) distribution as follows: , where the DE distribution places a higher density at zero and thicker tails, inducing stronger shrinkage of estimates for markers with relatively small effect and less shrinkage for markers with substantial effect.The residual variance (σ 2 e ) is specified as a scaled inverse chi-squared prior density, with degrees of freedom df e and scale parameter S e .
BayesA method 14,33 considers Gaussian distribution with null mean as prior for SNP marker effects (u w ) , and a SNP marker-specific variance ( σ 2 w ).The variance associated with each marker effect assumes a scaled inverse chi-square prior distribution, p σ 2 w = χ −2 σ 2 w |df u , S 2 u , with degrees of freedom ( df u ) and scale parameter ( S 2 u ) treated as known 14 .Thus, BayesA places a t-distribution for the marker's effects, i.e., p u w |df u , S 2 = t 0, df u , S 2 u , providing a thicker-tail distribution compared to the Gaussian, allowing a higher probability of moderate to large SNP effects.
BayesB assumes that a known proportion of SNP markers have a null effect (i.e., a point of mass at zero), and a subset of markers with a non-null effect that follow univariate t-distributions 3,12 , as follows: where π is the proportion of SNP markers with null effect, and 1 − π is the probability of SNP markers with non-null effect contributing to the variability of the FE-related trait 3 .Thus, the prior distribution assigned to SNP with non-null effects is a scaled inverse chi-square distribution.
BayesC method 34 assumes a spike-slab prior for marker effects, which refers to a mixture distribution comprising a fixed amount with probability π of SNP markers have a null effect, whereas a probability of 1 − π of markers have effects sampled from a normal distribution.The prior distribution is as follows: www.nature.com/scientificreports/Where σ 2 w is the common variance for marker effect, df u and S 2 B are the degrees of freedom and scale parameter, respectively, p 0 and π 0 [0,1] are the prior shape parameters of the beta distribution.

Machine learning techniques
Two machine learning (ML) algorithms were applied for genomic prediction: Multi-layer Neural Network (MLNN) and support vector regression (SVR).The ML approaches were used to alleviate the standard assumption adopted in the linear methods, which restrict to additive genetic effects of markers without considering more complex gene action modes.Thus, ML methods are expected to improve predictive accuracy for different target traits.To identify the best combination of hyperparameters (i.e., parameters that must be tuned to control the learning process to obtain a model with optimal performance) in the supervised ML algorithms (MLNN and SVR), we performed a random grid search by splitting the reference population from the forward scheme into five-folds 35 .

Multi-layer neural network (MLNN)
In MLNN, handling a large genomic dataset, such as 305,128 SNPs, is difficult due to the large number of parameters that need to be estimated, leading to a significant increase in computational demand 36 .Therefore, an SNP pre-selection strategy based on GWAS results in the training population using an MTGBLUP method (Fig. 1A) was used to reduce the number of markers to be considered as input on the MLNN.In addition, this strategy can remove noise information in the genomic data set.In this study, the traits displayed major regions explaining a large percentage of genetic variance, which makes using pre-selected markers useful 37 .
The MLNN model can be described as a two-step regression 38 .The MLNN approach consists of three different layer types: input layer, hidden layer, and output layer.The input layer receives the input data, i.e., SNP markers.The hidden layer contains mapping processing units, commonly called neurons, where each neuron in the hidden layer computes a non-linear function (activation) of the weighted sum of nodes on the previous layer.Finally, the output layer provides the outcomes of the MLNN.Our proposed MLNN architecture comprises two fully connected hidden layers schematically represented in Fig. 1B.The input layer in MLNN considered SNP markers that explained more than 0.125% of the genetic variance for FE-related traits (Fig. 1A; ~ 15 k for ADG and DMI, and ~ 16 k for FE and RFI).The input covariate X = {x p } contains pre-selected SNP markers (p) with a dimension Nxp (N individuals and p markers).The pre-selected SNP markers are combined with each k neuron (with k = 1, …, Nr) through the weight vector ( W ) in the hidden layer and then summed with a neuron-specific bias ( b k ) for computing the linear score for the neuron k as:Z [1] i = f (b k [1] + XW [1] ) (Fig. 1B).Subsequently, this linear score transformed using an activation function f (.) to map k neuron-specific scores and produce the first hidden layer output ( f z 1,i ).In the second-hidden layer, each neuron k receives a net input coming from hidden layer 1 as: www.nature.com/scientificreports/Z [2] i = b k [2] + Z [1] i W [2] , where W [2] represents the weight matrix of dimension k x k (k-number of neurons) connecting the Z [1] i into the second hidden layer, and b k [2] is a bias term in hidden layer 2.Then, the activation function is applied to map the kth hidden neuron unit in the second hidden layer and generate the output layer as V 2,i = f z 2,i .In the MLNN, a hyperbolic tangent activation function ( tanh(x) = e x − e −x /e x + e −x ) was adopted in the first and second layers, providing greater flexibility in the MLNN 39 .
The prediction of the adjusted FE-related trait was obtained as follows 38 : where y * represents the target adjusted feed efficiency-related trait for the ith animal; k the number of neurons considered in the model and assumed the same in the first and second layer; W 0 represents the weight from the k neuron in layer 2, b is related to the bias parameter.The optimal weights used in MLNN were obtained by minimizing the mean square error of prediction in the training subset 40 .
The MLNN model was implemented using the R package h2o (https:// github.com/ h2oai/ h2o-3), with the random grid search using the h2o.grid function (https:// cran.r-project.org/ web/ packa ges/ h2o) to determine the number of neurons to maximize the prediction accuracy.We used the training population split into fivefold to assess the best neural network architecture and then apply it in the disjoint validation set 41,42 .We considered a total of 1000 epochs 36 , numbers of neurons ranging from 50 to 2500 with intervals of 100, and applied a dropout ratio of 0.2 and regularization L 1 and L 2 parameters as 0.0015 and 0.0005, respectively.In this framework, the MLNN was performed using two hidden layers of neural networks with the number of neurons (k) of 750 for ADG, 1035 for DMI, 710 for FE, and 935 for RFI obtained during the training process.

Support vector regression
Support vector regression (SVR) is a kernel-based supervised learning technique used for regression analysis 43 .In the context of GS, the SVR uses linear models to implement nonlinear regression by mapping the predictor variables (i.e., SNP marker) in the feature space using different kernel functions (linear, polynomial, or radial basis function) to predict the target information, e.g., adjusted phenotype the GS 44 .SVR can map linear or nonlinear relationships between phenotypes and SNP markers depending on the kernel function.The best kernel function mapping genotype to phenotype (linear, polynomial, and radial basis) was determined using the training subset split into fivefold.The radial basis function (RBF) was chosen as it outperformed the linear and polynomial (degree equal 2) kernels in the training process, increasing 8.25% in predictive ability and showing the lowest MSE.
The general model for SVR using a RBF function can be described as 38,45 : y * i = b + h(m) T w + e, where h(m) T represents the kernel radial basis function used to transform the original predictor variables, i.e.SNP marker information ( m ), b denotes the model bias, and w represents the unknown regression weight vector.In the SVR, the learn function h(m) T was given by minimizing the loss function.The SVR was fitted using an epsilon-support vector regression that ignores residual absolute value ( y * i − y * i ) smaller than some constant (ε) and penalize larger residuals 46 .
The kernel RBF function considered in the SVR follows the form: h(m) T = exp −γ�m i − m j � 2 , where the γ is a gamma parameter to quantity the shapes of the kernel functions, m and m i are the vectors of predictor vari- ables for labels i and j.The main parameters in SVR are the cost parameter ( C ), gamma parameter ( γ ), and epsilon ( ǫ ).The parameters C and ǫ were defined using the training data set information as proposed by Cherkasky and Ma 47 : C = max y * + 3σ y * , y * − 3σ y * and ǫ = 3σ y * ln(n)/n , in which the y * and σ y * are the mean and the standard deviation of the adjusted FE-related traits on the training population, and n represents the number of animals in the training set.The gamma (γ) was determined through a random search of values varying from 0 to 5, using the training folder split into fivefold.The better-trained SVR model considered the γ parameter of 2.097 for ADG, 0.3847 for DMI, 0.225 for FE, and 1.075 for RFI.The SVR was implemented using the e1071 R package 48 .

Predictive ability
Prediction accuracy (acc) of the different statistical approaches was assessed by Pearson's correlation between adjusted phenotypes ( y * ) and their predicted values ( y * i ) on the validation set, and root mean squared error (RMSE).The prediction bias was assessed using the slope of the linear regression of y * i on y * , for each model.The Hotelling-Williams test 49 was used to assess the significance level of the difference in the predictive ability of Bayesian methods (BayesA, BayesB, BayesC, BL, and BRR), MTGBLUP, and machine learning (MLNN and SVR) against STGBLUP.The similarity between the predictive performance of the different models was assessed using Ward's hierarchical clustering method with an Euclidian distance analysis.The relative difference (RD) in the predictive ability was measured as RD = (r m −r STGBLUP ) r STGBLUP × 100 , where r m represents the acc of each alterna- tive approach (SVR, MLNN, and MTGBLUP, or Bayesian regression models-BayesA, BayesB, BayesC, BL, and BRR), and r STGBLUP is the predictive ability obtained using the STGBLUP method.

Ethical approval
The animal procedures and data sampling presented in this study were approved and performed following the Animal Care and Ethical Committee recommendations of the São Paulo State University (UNESP), School of Agricultural and Veterinary Science (protocol number 18.340/16).

Description of phenotypic and genotypic information
Principal component (PC) analysis identified some genetic stratification among the selection herds in the population, with the first two PCs explaining 3.27% of genetic variation (Additional File1: Figure S1).As expected, two selection herds in the IZ population (NeS and NeT) exhibited more genetic similarity, whereas the NeC herd was a completely separated cluster.The heritability estimates for FE-related traits were moderate, ranging from 0.21 ± 0.029 for FE to 0.40 ± 0.039 for DMI (Table 1).These heritability estimates indicated that a substantial component of the phenotypic variation in FE-related traits is due to additive genetic effects.The genetic correlations between FE-related traits were from low (r = 0.22 for RFI x ADG) to high moderate (r = 0.66 for RFI x DMI), and the FE showed a negative correlation with DMI and RFI (Supplementary Figure S2).

Predictive ability of GS approaches
Prediction accuracy comparison used GBLUP as the benchmark against ML methods (SVR and MLNN), MTG-BLUP, and Bayesian regression techniques (BayesA, BayesB, BayesC, BRR, and BL) (Table 2).The prediction accuracy for forward validation is shown in Table 2.The predictive correlation obtained for the various models varied from moderate to high, ranging from 0.53 (BayesA) to 0.67 (SVR) for ADG, 0.54 (BayesA) to 0.62 (MTG-BLUP and SVR) for DMI, 0.53 (STGBLUP) to 0.64 (MTGBLUP and SVR) for FE, and 0.62 (STGBLUP) to 0.69 (SVR) for RFI (Table 2).The accuracy obtained for the various models was directly proportional to the trait heritability (Supplementary Figure S3).Overall, the ML techniques (MLNN and SVR) and multi-trait GBLUP resulted in greater prediction accuracies against the standard approaches for GS, i.e., STGBLUP and Bayesian approaches (Table 2).
Table 1.Descriptive statistics for phenotype and adjust phenotype, estimates of additive genetic variances ( σ 2 a ), heritabilities ( h 2 ) with its respective standard error (SE) and the avarege accuracy of genomic breeding value ( Acc GEBV ) of all animals with standard deviation (SD) for feed efficiency-related traits, using the multi- trait GBLUP.ADG: average daily gain, FE: feed efficiency, RFI: residual feed intake, DMI: dry matter intake and Adjusted phenotype: phenotype information corrected for fixed effects considered in the model.Table 2. Accuracies of genomic prediction and root mean square error (RMSE) assessed for feed efficiencyrelated traits in testing animals using a forward validation scheme obtained with different parametric and machine learning methods. 1 ADG: average daily gain, DMI: dry matter intake, FE: feed efficiency, and RFI: residual feed intake, r: training-cor(y * , y * i ) in training population, Accuracy-c or(y * , y * i ). 2 STGBLUP: single trait GBLUP, BayesA: Bayesian A, BayesB: Bayesian B, BayesC: Bayesian C, BL: Bayesian Lasso, BRR: Bayesian ridge regression, MTGBLUP: multi-trait GBLUP, MLNN: Multi-layer neural networks, and SVR: support vector machine regression using a radial basis kernel.Statistically significant differences between the predictive ability of each method across the herds compared to the standard model (STGBLUP) were *p-value < 0.05, ** p-value < 0.01 and ***p-value < 0.005.

Bias and residual parameters of the GS approach
The slope coefficient of regression of the predicted and observed (adjusted) phenotype was used to assess the extent of prediction bias (Fig. 3A).Unbiased GS approaches are expected to show a regression slope near 1.The results show that the slope coefficients differed significantly from 1 for the Bayesian methods (BayesA, BayesB, BayesC, BRR, and BL) and STGBLUP (Fig. 3A), so the predictions were empirically biased for these models.On the other hand, slope coefficients for SVR and MTGBLUP differed only slightly from 1, so the predictions were empirically unbiased for these approaches, while MLNN showed biased predictions only for ADG (Fig. 3A).Using a hierarchical cluster based on accuracy and slope regressions across traits, we observed that alternative approaches (MTGBLUP, MLNN, and SVR) were similar to each other but different from Bayesian regression methods and STGBLUP (Fig. 3B).However, while Bayesian models and STGBLUP tended to be similar, while Bayesian Lasso (BL) showed a greater distance from the others due to its lower predictive ability and highest slope regressions (Fig. 3B).
Assessment of prediction quality through RMSE indicated that ML techniques and MTGBLUP led to considerable reductions in the predictive error compared to STGBLUP (Table 2).Concerning RMSE, the alternative approaches (MTGBLUP, MLNN, and SVR) had the best performance for all traits, with a higher reduction in RMSE compared to STGBLUP and Bayesian regression, varying from 21 to 29% with MTGBLUP, 10% to 14% with MLNN, and 21% to 30% with SVR (Table 2).On the other hand, BL had the highest RMSE value compared with the other approaches, which might be caused by a more biased prediction (Table 2).Comparing the residual values by the probability distribution indicated that SVR and MTGBLUP had a higher probability of lower error than any other approach (Fig. 4).The BL approach showed the highest probability of large residual values, suggesting more biased predictions (Fig. 4).

Discussion
Genomic prediction of complex traits has been based mainly on parametric methods, which generally consider only additive genetic effects.However, recently, alternative models such as ML or MTGBLUP have been considered as promising approaches for GS due to the potential to capture nonadditive SNP effects for ML or to leverage the genetic correlation between traits for MTGBLUP, and hence improving prediction accuracy 4,18,50,51 .This study compared the prediction accuracy of commonly applied parametric methods (STGBLUP and Bayes regression) with MTGBLUP and nonparametric ML methods (MLNN and SVR) for FE-related traits.We used adjusted phenotypes as the target variables to reduce any noise from contemporary group effects in the phenotype while considering all gene actions.Our results show that the nonlinear prediction methods (MLNN and SVR) Figure 2. Relative difference in predictive ability using Bayesian methods (BayesA-Bayesian A, BayesB-Bayesian B, BayesC-Bayesian C, BL-Bayesian Lasso and BRR-Bayesian ridge regression), multi-trait GBLUP (MTGBLUP), and machine learning (MLNN-Multi-layer neural networks, and SVR-support vector machine regression) against to single-trait GBLUP (STGBLUP).and MTGBLUP outperform the other methods for FE-related traits in Nellore cattle (Table 2).However, when the GEBV is used as the target variable in a breeding context, models accounting for exclusively additive effects attributed to markers generally perform better in traits with purely additive inheritance 15 .Nevertheless, Liang et al. 52 found that combining three machine learning approaches using a stacking ensemble learning framework can accurately predict GEBV for different traits with varying genetic architecture in beef and dairy cattle.Huang and Mackay 53 have suggested that the estimated additive variance resulting from the similarity between relatives may include a portion of the nonadditive variance.This could explain why nonlinear methods perform better, even if the target variable is expected to depend only on additive effects (GEBV).Depending on the complexity of the trait and the prediction scenario, more sophisticated models help to increase the prediction accuracy.In this situation, using simple models may miss essential data features and cause losses of prediction accuracy for complex traits, whereas nonlinear models usually give better prediction accuracy 19 .
Previous studies reported prediction accuracies of GS for FE-related traits in Nellore cattle ranging from 0.06 to 0.54, obtained with single-step GBLUP and BayesCπ 6 .Our results indicated moderate prediction accuracy and higher than those obtained in other studies for FE-related traits 6,54,55 .In this study, an important factor leading to increased prediction accuracy is the relationship between the reference and testing population (Supplementary Figure S4), which has been pointed out as an important factor to enhance prediction accuracy 5,6,[55][56][57] .
The usefulness of the statistical methods for GS depends on the extent to which these approaches are able to cope with the genetic architecture underlying the target trait and handle the problems associated with SNPs in high linkage disequilibrium 38 .Performing GS using models with less strong assumptions regarding the SNP effects (MLNN and SVR) or considering the correlation structures between the traits (MTGBLUP) outperformed the other models (STGBLUP and Bayesian regression).It has been reported that ML approaches lead to more accurate predictions when the target traits are affected by nonadditive effects due to their ability to account for these complex effects without explicit parametric representation in the models 58 .
The MLNN captures interactions between SNP markers and allows weight differentially in each neuron, reflecting its importance for output prediction 17,59 .On the other hand, the SVR approach uses a RBF to assign flexibility to markers, attempting to match correctly their contributions to the variation of the target trait 60,61 .www.nature.com/scientificreports/Although both the SVR and MLNN methods potentially capture nonadditive relationships between phenotypes and SNP markers, SVR using RBF shows advantages over MLNN in terms of more flexibility to handle a smaller number of animals and a larger number of genetic markers; MLNN requires a larger training population to learn from the input data satisfactorily provide accurate predictions with complex genetic architecture of traits 17,62 .
Another factor that could explain the slight differences in results among the approaches is that FE-related traits may be indeed affected by several QTLs with minor effects distributed across the genome 63,64 .This may have favored the models based on the infinitesimal assumption (BRR and STGBLUP) and limited the performance of the Bayesian models, leading to roughly equivalent prediction accuracies.The usefulness of the Bayesian regression approaches, which perform variable selection and/or shrinkage of markers, occurs when few QTLs with major effects explain a large proportion of the trait genetic variance 10,65,66 .Several studies have also reported that these traditional parametric models, with conceptual differences in prior assumptions on marker effects, showed similar or slight differences in predictive ability compared to STGBLUP 4,[67][68][69] .
Statistical methods capturing additive and epistatic effects have recently been reported to perform better than a simple additive model for the genomic prediction of traits with complex inheritance 70,71 .In order to take advantage of nonlinear relations between genomic information and phenotypes, nonparametric methods as machine learning techniques without prior knowledge about the distribution of the SNP marker effects have been exploited 18,20,51,72 .These methods can potentially capture nonadditive marker effects, improving the prediction quality compared to parametric methods 73 .Our results showed that machine learning techniques (MLNN and SVR) and MTGBLUP outperformed conventional methods (STGBLUP and Bayesian models) in predictive ability with unbiased predictions for the FE-related traits (Table 2 and Fig. 3), potentially capturing nonadditive genetic effects.Other studies have also shown that SVR, MLNN, and multi-trait GBLUP are promising approaches capable of outperforming traditional parametric methods commonly used for complex traits controlled by many genes with possible epistatic and pleiotropic effects 4,20,50,[74][75][76] .
The performance of the MLNN approach significantly improved the prediction accuracy by 9.8% and reduced the RMSE by 11% on average, compared to STGBLUP (Table 2 and Fig. 2), indicating its potential usefulness for genome-enabled prediction of FE-related traits in Nellore cattle.These results agree with other studies that observed a reasonable increment in prediction accuracy, although there was no consistent evidence indicating that MLNN can outperform linear models 71,[77][78][79] .The relative performance of MLNN vs. traditional methods depends on the association between the complexity of the genetic architecture of the target trait and the MLNN architecture.Bellot et al. 80 observed that deep learning techniques such as Multi-layer Perceptrons and Convolutional Neural Networks exhibited a slight reduction in predictive ability compared to Bayesian regression methods (BayesB and BRR) for traits with an additive and polygenic genetic background in humans.
Based on the overall predictive ability, the MTGBLUP and SVR methods outperformed all the other methods (MLNN, S-TGBLUP, and Bayesian regressions).A favorable attribute of SVR algorithms for GS is related to their ability to detect the most informative SNP markers associated with the target phenotypes and use these

Figure 1 .
Figure 1.(A) Manhattan plot for percentage of genetic variance explained by SNP-marker estimated through multi-trait GWAS in training population to be used as pre-selection strategies for multi-layer neural network.(B) General representation of neural networks with two hidden layers used to model nonlinear dependencies between trait and SNP marker information.The input layer ( X = x i,p ) considered in the neural network refers to the SNP marker information (coded as 0, 1, and 2) of the ith animal.The selected node represents the initial weight ( W = w p ), assigned as random values between -0.5 and 0.5, connecting each input node to the first hidden layer and in the second layer the w up refers to the output weight from the first hidden layer, b represents the bias which helps to control the values in the activation function.The output ( y ) layer represents a weighted sum of the input features mapped in the second layer.

Figure 3 .
Figure 3. Barplot of the bias of prediction determined as the slope of the linear regression of predicted adjusted phenotype on adjusted phenotype for feed efficiency-related trait using different genomic prediction approaches (A) and Ward's hierarchical clustering of models based on predictive ability and bias for each model (B).

Figure 4 .
Figure 4. Residual diagnostics comparison for forward validation assessed by the absolute residual value for machine learning (SVR and MLNN) and parametric approaches (BayesA, BayesB, BayesC, BL, BRR, STGBLUP, and MTGBLUP) the genomic prediction approach for the feed efficiency-related traits (ADG, DMI, FE, and RFI) in Nellore cattle.